# ifndef EQUATIONSOLVER_H
# define EQUATIONSOLVER_H

# define _USE_MATH_DEFINES
# include<cmath>
# include<limits>

const double _eps = 50 * std::numeric_limits<double>::epsilon();

/**
* @brief 定义一个抽象函数类，用于规范仿函数的接口以及给派生类提供求导数运算。
* 
*/
class Function
{
public:
    virtual double operator()(double _x) = 0;

    // 求导数运算
    virtual double diff(double _x)
    {
        return ((*this)(_x + _eps) - (*this)(_x)) / _eps;
    }
};


/**
* @brief 定义一个抽象类规范求解器的求解函数接口
*
*/
class EquationSolver
{
public:
    virtual double solve() = 0;  // 虚类，不得实体化
};

/**
* @brief 实现二分法(Bisection method)求零点的求解器
*
*/
class Bisection : public EquationSolver
{
private:
    double a, b, delta, eps;
    int M;
    Function & func;

public:
    Bisection(double _a, double _b, double _delta, double _eps, 
              int _M, Function & _f): func(_f)
    {
        a = _a; b = _b; delta = _delta;
        eps = _eps; M = _M;
    };
    ~Bisection(){};

    // 求解函数
    double solve()
    {
        double h = b - a;
        double u = func(a);
        double c = 0, w = 0;
        for (int k = 1; k <= M; k++)
        {
            h = h/2;
            c = a + h;
            w = func(c);
            if (fabs(h) < delta || fabs(w) < eps)
                break;
            else if ((w>0 && u>0) || (w<0 && u<0))
                a = c;
        }
        return c;
    };
};

/**
* @brief 实现牛顿迭代法(Newton's method)求零点的求解器
*
*/
class Newton : public EquationSolver
{
private:
    double x0, eps;
    int M;
    Function & func;

public:
    Newton(double _x0, double _eps, int _M, Function & _f): func(_f)
    {
        x0 = _x0; eps = _eps; M = _M;
    }
    ~Newton(){};

    // 求解函数
    double solve()
    {
        double x = x0, u;
        for (int k=0; k<M; k++)
        {
            u = func(x);
            if (fabs(u)<eps) break;
            x = x - u/func.diff(x);
        }
        return x;
    }
};

/**
* @brief 实现割线法(Secant method)求零点的求解器
*
*/
class Secant : public EquationSolver
{
private:
    double x0, x1, delta, eps;
    int M;
    Function & func;

public:
    Secant(double _x0, double _x1, double _delta, double _eps,
           int _M, Function & _f): func(_f)
    {
        x0 = _x0; x1 = _x1; delta = _delta;
        eps = _eps; M = _M;
    }
    ~Secant(){};

    // 求解函数
    double solve()
    {
        double xn = x1, xn1 = x0, s;
        double u = func(xn);
        double v = func(xn1);
        for (int k=2; k<=M; k++)
        {
            if(fabs(u)>fabs(v))
            {
                double t = xn; xn = xn1; xn1 = t;
                t = u; u = v; v = t;
            }
            s = (xn-xn1)/(u-v);
            xn1 = xn;
            v = u;
            xn = xn - u*s;
            u = func(xn);
            if (fabs(xn-xn1)<delta || fabs(u)<eps)
                break;
        }
        return xn;
    }
};

# else
// Do nothing.
# endif